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POINT  RESPONSE  FUNCTIONS  FOR  UNFILLED  OPTICAL  ARRAYS 


ABSTRACT 

Two  problems  relating  to  the  design  of  arrays  for  Imaging  complex  fields 
are  the  sidelobe  levels  and  the  energy  within  those  sidelobes.  While  design¬ 
ers  often  know  the  performance  of  a  filled  aperture,  the  much  larger  number  of 
parameters  involved  in  array  design  lead  to  great  difficulty  evaluating  the 
performance  level.  This  report  describes  a  family  of  optical  arrays  whose 
parameters  were  evaluated  by  a  computer  analysis  system.  For  simplicity  all 
the  arrays  considered  were  in  the  same  family  (a  circular  ring  of  evenly 
spaced  array  elements) ,  so  that  comparisons  could  be  made  as  the  number  of 
array  elements  was  varied.  Arrays  with  3,  7,  15,  and  27  elements  were  con¬ 
sidered  in  detail,  including  the  proper  corrections  for  both  finite  element 
size  and  finite  bandpass. 

We  conclude  that  a  spaceborne,  unfilled,  minimum-redundancy  array  can  be 
built  for  earth  viewing  that  has  performance  similar  to  a  filled  aperture 
significantly  reducing  the  cost  of  space  based  imaging  systems.  The  sidelobe 
levels  and  energy  in  the  sidelobes  can  be  controlled  through  the  appropriate 
selection  of  design  parameters.  However  a  final,  detailed  design  will  depend 
upon  target  characteristics  and  cost-benefit  tradeoffs,  as  well  as  upon  the 
center  wavelength,  bandpass,  desired  resolution,  and  the  allowable  sidelobe 
levels.  The  analysis  methods  we  have  developed  have  general  applicability  to 
the  design  problem. 


INTRODUCTION 

The  use  of  interferometer  systems  in  applications  where  they  previously 
could  not  be  used  is  now  possible  because  of  great  advances  in  both  inter¬ 
ferometer  and  image  processing  techniques  in  the  past  decade.  This  is  par¬ 
tially  due  to  better  signal  processing  hardware  and  software  and  because 
imaging  problems  arose  which  needed  novel  solutions.  While  interferometers 
provide  increased  resolution  at  low  cost,  their  increased  sidelobe  levels  can 
create  major  problems,  especially  for  an  imaging  system.  In  a  simple  field 
where  the  received  signal  is  localized  to  a  small  percentage  of  the  total 
image,  high  sidelobes  can  be  tolerated  and  elementary  methods  restore  the 
original  brightness  distribution  satisfactorily.  However,  in  the  case  of  a 
down-looking  system  from  high  altitudes  or  space,  the  image  brightness  is 
complicated;  there  is  information  in  every  part  of  the  image,  much  of  which  is 
irrelevant  and  mu3t  be  removed  to  recognize  a  target.  Removing  the  background 
is  almost  impossible  if  the  point  response  function  of  the  system  has  high 
sidelobes  extending  across  the  image,  because  objects  are  not  localized.  We 
describe  an  optical  array  that  has  low  sidelobe  levels  and  a  high  percentage 
of  the  power  in  the  main  beam. 

Of  possible  greater  importance,  this  technical  memorandum  also  describes 
a  new  computer  analysis  system  for  examining  the  sidelobe  levels  and  their 
effect  upon  an  image.  A  family  of  spaceborne  optical  arrays  have  been  ex¬ 
amined  to  demonstrate  these  techniques.  It  should  be  emphasized  that  there  is 
nothing  special  in  the  demonstration  array.  Using  this  analysis  system,  a 
completely  different  configuration  could  be  thoroughly  examined  with  a  few 
weeks  of  work. 

Manuscript  approved  March  12,  1984. 
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Radio  astronomers  are  not  able  to  build  filled  apertures  to  achieve  the 
resolution  they  desire.  Since  adequate  sensitivity  can  be  achieved  with  an 
unfilled  aperture,  they  have  built  interferometer  systems  (ref.  1,2).  They 
have  had  the  problem  of  high  sidelobes  for  many  years  because  of  practical 
limits  to  the  number  of  elements  in  any  array.  As  a  result,  extensive  soft¬ 
ware  has  been  developed  for  calculating  point  response  functions  for  realistic 
cases  using  arbitrarily  located  elements.  The  basic  approach  of  this  memo¬ 
randum  is  therefore  to  apply  some  of  these  special  computer  programs  to  the 
specific  case  of  an  optical  array,  writing  new  programs  to  aid  the  interfacing 
or  to  do  peculiar  cases. 

We  first  review  the  optical  theory  which  leads  to  the  point  response 
function  of  addition  correlators.  This  is  a  different  situation  from  the  more 
familiar  multiplication  correlators  used  at  radio  wavelengths.  Before  exami¬ 
ning  complex  interferometers,  the  radial  distribution  of  the  point  response 
function  for  a  filled  aperture  is  derived  and  compared  to  a  circularly  sym¬ 
metric  interferometer.  Finally  the  discrete  element  broad  bandpass  case  is 
synthesized. 


THE  THEORY  OF  ADDITION  CORRELATORS 

The  first  step  in  understanding  optical  interferometers  is  to  review  the 
physics  of  addition  correlators.  At  the  present  it  is  not  technically  poss¬ 
ible  to  make  multiplying  optical  interferometers  similar  to  conventional  radio 
interferometers.  Following  the  basic  development  of  Born  and  Wolf  (ref.  3), 
we  see  that  for  a  plane  wave  the  total  intensity  is  related  to  the  electric 
vector  by: 


I  = 


(1) 


—  i  cot 

where  the  electric  vector,  E(r,t)  =  Re  {A(r)e  } ,  is  a  function  of  ^pace  + 
and  time.  Now  consider  the  superposition  of  two  monochromatic  waves  E^  and  E 2 
at  a  point  in  space,  P 

E  =  E1  +  E2 


or 


<fa>  =  <f1*>  +  <§2*>  +  2  <E1  *  S2> 


(2) 


Hence  the  total  intensity,  I,  at  point  P  is 


1=1. 


J2  +  J12 


where  we  have  made  use  of  the  definitions 


(3) 


I1  =  <E1l>  ;  I2  =  <E2*>;  J12  =  2  <E,  •  E2> 

Notice  that  the  intensities  I  ,  I_,  and  I  are  positive  definite  by  this  for¬ 
mulation,  but  no  such  requirement^exists  for  J Radio  astronomers  are  able 
to  measure  the  interference  terra  J12  separately^  eq.  3  because  they  use  a 
multiplication  correlator  instead  or  the  addition  correlator  as  was  used  in 
obtaining  eq.  3*  As  a  result  radio  astronomers  can  obtain  negative  sidelobes 
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and  negative  portions  of  their  images,  while  that  is  impossible  for  the  addi¬ 
tion  correlator. 

Expanding  this  analysis  to  n  elements  is  relatively  easy. 

I  =  <E  *>  +  <1  *>  ♦  ...  +  <E  *> 

12  n 

+  2  <E  *  E„>  +  ...  +  2  <E  *  E  > 

12  In 

+  ...  +  2  <E  .  *  E  >  (4) 

n— I  n 

The  sun  of  all  the  cross  terms  J  divided  by  the  number  of  cross  terms  is 
what  radio  astronomers  call  the  ^Synthesized  beam".  There  are  standard  com¬ 
puter  programs  for  computing  the  synthesized  beam,  so  it  is  easy  to  analyze 
the  characteristics  of  an  optical  array  using  the  synthesized  beam  as  a  start¬ 
ing  point.  For  an  interferometer  with  n  identical  elements  there  are  n(n-1)/2 
cross  terms,  and  eq.  4  simplifies  to: 

I  s  n  (I|  +  (n-1)  [Synthesized  Beam])  (5) 

This  result  is  for  infinitesimally  small  point  interferometer  elements  and 
monochromatic  waves.  For  finite  elements  the  equation  must  be  multiplied  by 
the  individual  element  response  which  is  an  envelope  or  gain  function  for  the 
intensity  as  given  in  eq.  5.  For  a  finite  bandpass  the  fringe  functions  add 
out  of  phase  so  as  to  apply  a  taper  to  the  synthesized  beam  for  angles  off 
axis  (See  Appendix  A).  Therefore  a  bandpass  effect  is  applied  only  to  the 
second  term  on  the  right  hand  side  of  eq.  5.  If  we  denote  the  (normalized) 
response  of  the  individual  elements  by  A(x,y)  and  the  bandpass  effect  by 
B(x,y)  we  obtain, 

I(x,y)  =  n  A(x,y)  [I^x.y)  +  (n-1)  B(x,y)  b(x,y)]  (6) 

where  we  have  finally  denoted  the  synthesized  beam  by  b(x,y).  We  will  be  able 
to  examine  these  terms  more  closely  by  using  several  special  cases.  A  most 
important  special  case  is  when  the  number  of  elements  approaches  infinity  and 
the  aperture  is  filled.  This  limiting  case  is  examined  next. 


THE  POINT  RESPONSE  OF  A  FILLED  APERTURE 

Interferometers  do  not  measure  the  brightness  of  the  image  directly, 
but  measure  the  Fourier  transform  of  the  image  which  is  the  image  visibility 
V(u,v).  In  the  case  of  a  point  source  at  the  (x,y)  origin,  the  visibility  is 
a  constant  vector  with  zero  phase  and  an  amplitude  equal  to  the  flux  density 
of  the  point  source.  The  synthesized  bean  of  an  array  of  interferometers  can 
be  written  as 


b(x,y)  [*  *]  S(u,v)  w(u,v)  (7) 

where 

b(x,y)  =  synthesized  beam. 

S(u,v)  =  sampling  function  that  represents  the  distribution  of 
the  data  points  in  the  (u,v)  plane. 
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w(u,v)  s  weighting  function  of  the  data  including  any  applied 
tapering.  S(u,v)  w(u,v)  is  the  optical  modulation 
transfer  function. 

[*  *]  denotes  Fourier  transformation. 

The  function  S(u,v)w(u,v)  is  entirely  real  so  b(x,y)  is  symmetrical. 

In  real  situations  there  is  finite  bandwidth.  The  mathematics  for  two 
element  Interferometers  with  finite  bandpass  and  for  a  simple  5  percent  rec¬ 
tangular  bandpass  array  are  derived  in  Appendix  A.  The  two  element  result  is 
a  reduction  in  the  fringe  amplitude  (reduction  in  the  sidelobe  level)  as  the 
distance  from  the  delay  center  increases  with  the  functional  relationship  of 
sin(x)/x  where  this  x  is  *  times  angular  distance  divided  by  baseline.  The 
general  case  of  discrete  arrays  with  large  bandpass  is  not  solved  in  closed 
form;  for  the  calculations  of  this  report  we  sum  the  n(n-1)/2  sin(x)/x  terms 
directly  using  a  computer. 

The  monochromatic  point  response  function  of  a  circular  filled  aperture 
is  then  simply  the  Fourier  Transform  of  a  uniform  disk,  the  first  order  Bessel 
function  J.(vr/a)/(r/a)  where  a  is  the  diameter  of  the  aperture.  This  radial 
dependency1 is  shown  in  Fig.  1  where  we  have  assumed  a  diameter  of  10  meters 
and  a  center  wavelength  of  500  nm.  When  the  effects  of  the  finite  rectangular 
bandpass  from  350  nm  to  650  nm  are  taken  into  account,  the  point  response 
function  is  effectively  multiplied  by  the  delay  response,  shown  in  Fig.  2a. 

The  broadband  result  is  shown  in  Fig.  2b. 


Figure  1.  Point  response  function  of  a  monochromatic  filled  aperture  5  meters 
in  radius. 
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THE  CIRCULARLY  SYMMETRIC  INTERFEROMETER 


We  next  calculate  the  results  for  the  case  of  an  infinite  number  of 
point  elements  arranged  in  a  circle.  This  important  case  is  solvable  in 
closed  form  (ref.  4)  and  the  synthesized  beam  b(x,y)  is  given  by 

b(x,y)  =  JQ(vr/a)  (8) 

where  a  is  used  to  denote  the  radius  of  the  circle.  This  function  (shown  in 
Fig.  3)  does  not  decrease  with  r  nearly  as  rapidly  as  J.(x)/x  which  we  found 
for  a  filled  aperture.  However,  in  realistic  cases  finite  element  size  and 
bandpass  provide  a  taper  to  lower  the  sidelobe  level  more  rapidly  than  eq.  8. 

Combining  eqs.  6  and  8  leads  to 

I(x,y)  =  A(x,y)  B(x,y)  JQ(*r/a)  (9) 

because  as  n  approaches  infinity  the  first  term  of  eq.  6  can  be  neglected. 

This  is  the  same  result  as  for  a  multiplication  correlator. 

At  a  distance  of  1  arcsecond  the  sidelobes  have  dropped  to  30  dB  down 

from  the  central  maximum,  but  are  50  dB  higher  than  for  the  circular  aperture 

case.  Because  these  profiles  are  well  known  and  understood,  it  is  easy  to  U3e 

spatial  filtering  to  make  the  ring  response  match  the  filled  aperture.  The 

matched  filter  is  [x  J  (x)/J,(x)]*. 

o  1 

If  the  effects  of  finite  element  size  are  included  into  the  profile 
shown  in  Fig.  3a,  the  function  drops  more  rapidly.  The  10  cm  elements  assumed 
here  have  a  minimum  at  a  full  width  half  maximum  of  1  arcsecond.  If  a  single 
detector  is  used  at  each  element  (as  in  a  radio  interferometer),  then  the 
finite  bandpass  response  will  severely  narrow  the  array  response,  providing  an 
additional  sidelobe  suppression.  The  combined  effects  of  finite  elements  and 
bandpass  are  shown  in  Fig.  3*>. 

THE  DISCRETE  ELEMENT  ANALYSIS 

The  problem  of  simulating  the  point  response  function  at  the  origin  of  a 
complex  array  with  finite  elements  and  a  broad  bandpass  is  the  special  case  of 
simulating  the  response  to  any  image.  To  allow  for  more  general  cases,  the 
best  approach  is  to  simulate  the  array  (using  software)  for  a  specific  array 
geometry,  wavelength,  bandpass,  etc.  We  chose  a  ring  geometry  for  the  ele¬ 
ments,  since  that  maximizes  the  two  dimensional  resolution  and  yields  a  rea¬ 
sonable  compromise  for  the  minimum  element  spacing.  Thus,  a  ring  is  a  system 
with  a  wide  range  of  spatial  frequencies.  We  analyzed  the  realistic  case  of 
finite  circular  elements  arranged  on  a  ring.  We  have  considered  in  some 
detail  ring  arrays  of  3,  4,  7,  9,  15,  27,  45,  and  99  elements,  but  we  have  not 
applied  the  full  details  of  the  analysis  in  all  cases.  As  the  number  of 
elements  in  the  ring  increases,  they  are  closer  spaced,  which  pushes  the  first 
grating  ring  further  out.  Within  this  ring  the  discrete  element  case  is 
closely  approximated  by  the  infinitesimal  ring  discussed  above. 

For  our  analysis,  we  have  assumed  a  geometry  consisting  of  a  10  meter 
array  operating  at  a  wavelength  of  500  nanometers.  The  analysis  includes  both 
point  elements  and  the  more  realistic  case  of  10  cm  elements.  For  broadband 
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Figure  3.  Point  response  function  for  the  circularly  symmetric  unfilled 
aperture,  (a) Monochromatic  case,  (b) Finite  bandpass  case. 


effects  we  assume  a  uniform  wavelength  response  between  350  and  650  nano¬ 
meters.  Unless  specifically  stated  however,  results  are  for  monochromatic 
light  at  500  nanometers.  The  analysis  assumes  the  beams  from  array  elements 
are  combined  coherently.  This  assumption  implies  that  the  interferometer  is 
built  away  from  the  degrading  effects  of  the  atmosphere.  System  noise  and 
other  operational  degradations  are  not  a  part  of  the  point  response  function, 
and  were  not  considered.  Conceptually,  the  interferometer  we  are  considering 
consists  of  an  array  of  elements  on  a  ring,  whose  individual  beams  are  brought 
to  a  common  focus. 

The  analysis  system  we  have  developed  is  a  collection  of  specialized 
programs  from  two  different  image  analysis  systems  with  the  addition  of  a  few 
programs  specially  written  for  this  project.  The  basic  flow  shown  in  Fig.  4 
is  to  assume  an  interferometer  configuration,  generate  fake  data  for  it  using 
the  Caltech  Very  Long  Baseline  Interferometer  (VLBI)  package  (ref.  5),  convert 
the  fake  data  to  the  NRAO  image  processing  package  AIPS  (ref.  6),  and  invert 
the  data  to  obtain  the  synthesized  beam.  This  is  the  synthesized  beam  for 
monochromatic  point  elements.  Using  the  interactive  data  language  IDL  (ref. 

7) ,  we  generate  the  element  responses  for  the  interferometer  and  convert  it 
into  AIPS  format  where  it  can  be  combined  with  the  synthesized  beam  and  plott¬ 
ed  in  a  variety  of  ways.  The  bandpass  correction  function  was  introduced 
separately  from  a  special-purpose  Fortran  program. 

Because  both  of  the  major  packages  we  use  are  organically  growing  at  the 
present  time,  a  detailed  step-by-step  description  of  the  method  would  be  obso¬ 
lete  when  the  next  release  is  received  within  the  month.  Of  more  value  is  a 
flow  description  which  leads  to  an  understanding  for  the  need  for  each  step. 
The  exact  modeling  done  here  is  necessarily  quite  complicated  and  requires 
many  stages  through  the  software. 

The  first  step  in  the  analysis  (see  Fig.  4)  was  to  create  an  artificial 
interferometer  data  set,  using  the  Caltech  VLBI  program  FAKE.  In  appendix  B 
are  the  3  files  used  to  run  this  program  in  the  case  of  the  15-element  array. 
The  parameters  in  the  FAKE.COM  file  were  selected  to  give  a  "noise-free"  data 
set  symmetric  about  the  V  axis  in  the  UV  plane.  The  file  STATIONS.DAT  con¬ 
tained  the  element  positions;  POINT. MOD  contained  the  parameters  for  a  dummy 
input  source  model  (for  the  purposes  of  our  analysis,  the  type  of  input  source 
model  used  was  irrelevant).  The  original  maps  made  from  the  model  (using  the 
data  generated  by  FAKE)  did  not  include  either  bandpass  or  finite  element 
effects  and  thus  could  not  be  used  directly.  The  model  source  data  is  dis¬ 
carded  following  the  UVSRT  program  step.  The  elements  were  located  on  a 
circle;  for  convenience  a  circle  of  latitude  was  used  with  the  frequency  being 
scaled  so  a3  to  give  a  circular  interferometer  with  the  desired  size  of  20 
million  wavelengths.  Further  explanation  of  the  parameters  used  for  FAKE  can 
be  found  in  the  documentation  for  the  Caltech  VLBI  package  (ref.  5). 

The  fake  data  was  then  translated  into  the  AIPS  format  using  the  AIPS 
program  T0AIP,  (or  FRMVB  for  more  than  27  elements),  sorted  as  necessary  using 
UVMAP ,  and  then  Fourier  transformed  by  UVMAP.  UVMAP  has  2  output  images:  the 
normalized  synthesized  beam  b(x,y)  for  the  array  under  consideration  and  the 
so-called  "dirty-map",  the  latter  of  which  was  discarded.  The  synthesized 
beam  is  the  point  response  function  for  the  monochromatic,  point-like  element 
case,  with  a  zero-level  offset  (see  eq.  6).  As  an  example,  the  synthesized 
beam  b(x,y)  for  a  15-element  array  is  presented  in  Fig.  5. 
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Figure  4.  Computer  analysis  logical  flow.  The  flow  chart  shows  the  relation¬ 
ship  between  the  various  programs  to  produce  point  response  functions  for 
finite  bandpasses  and  finite  elements.  The  numbered  points  on  the  flow  chart 
refer  to  intermediate  images  as  follows: 

(1)  FAKE.IBEM  is  the  monochromatic  synthesized  beam  b(x,y)  in 
equation  6. 

(2)  BW.MAP  is  the  bandpass  correction  function  B(x,y). 

(3)  PRF1.MAP  is  the  term  [I1(x,y)  +  n-1  B(xfy)  b(x,y)]  in 
equation  6,  which  represents  the  normalized  point  response 
function  for  an  array  of  point  elements  with  finite 
bandpass . 

(4)  MI. MAP  is  the  term  A(x,y)  in  equation  6,  representing  the 
beam  from  the  individual  elements  in  the  array. 

(5)  PRF.MAP  is  I(x,y)f  the  normalized  point  response 
function  for  finite  bandpass  and  finite  element  size. 
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Figure  5.  Synthesized  beam  function  for  the  discrete  element  unfilled  apcrt 
(monochromatic  case  and  point  elements').  The  array  consisted  o’"  1b  ’lament, 
arranged  on  a  ring.  The  scale  across  the  bottom  defines  color  as  i  functLc* 
brightness.  Thi3  map  corresponds  to  F^KE.IdEd  in  the  flow  chart,  of  trig. 
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The  taper  due  to  the  finite  bandpass  must  next  be  accounted  for.  For  the 
special  case  of  an  odd  number  of  elements  arranged  symmetrically  on  a  ring, 
program  BAND  (Appendix  B)  generates  the  bandpass  array  B(x,y)  by  which  the 
synthesized  beam  must  be  multiplied.  For  a  generalized  array  with  arbitrary 
element  locations,  BAND  could  not  be  used  as  it  stands,  but  would  have  to  be 
modified  to  read  the  UV  data  to  get  the  position  angle  and  scaling  for  the 
sin(x)/x  calculations.  An  example  of  the  B(x,y)  function  for  the  15-eleoent 
array  is  shown  in  Fig.  6. 

The  array  from  BAND  is  copied  into  an  existing  dummy  AIPS  image  file 
using  IDL  for  convenience  in  getting  the  correct  format.  To  use  IDL  in  this 
way,  one  first  opens  and  reads  the  data  file  into  an  internal  array  in  IDL. 
Next,  the  AIPS  dummy  file  is  opened  using  an  OPENU  statement  (see  the  IDL 
documentation,  ref.  7,  for  details).  Next,  a  variable  is  associated  with  the 
AIPS  file.  Then  the  internal  array  is  copied  into  the  associated  variable 
with  a  statement  like  A(0)=B.  Finally  the  AIPS  file  is  closed  in  the  standard 
way.  Care  must  be  taken  not  to  corrupt  the  normalization  of  the  AIPS  file  too 
badly,  or  unpredictable  things  can  happen.  Also,  IDL  must  write  the  file  with 
the  correct  data  type  needed  by  AIPS. 

The  multiplication  of  the  synthesized  beam  b(x,y)  by  the  bandpass  func¬ 
tion  B(x,y)  is  performed  by  the  AIPS  task  COMB.  In  the  same  operation,  the 
(n-1)  factor  for  the  product  (see  eq.  6)  is  applied  and  the  constant  I  term 
in  equation  (6)  is  added  to  produce  the  intermediate  map  PRF1 .  PRF1  is  the 
normalized  point  response  function  for  an  array  of  point  elements  having  a 
finite  bandpass.  PRF1  MAP  for  the  15-element  array  is  presented  in  Fig.  7. 
There  is  a  dramatic  improvement  in  the  sidelobe  level  due  to  B(x,y). 

The  final  step  in  producing  the  point  response  function  is  to  apply  the 
finite  element  correction  A(x,y).  This  correction  was  generated  using  the  IDL 
procedure  MIRROR  (see  appendix  B)  and  copied  into  AIPS  using  IDL  (see  above). 
COMB  was  then  used  to  multiply  and  renormalize  the  two  arrays  resulting  in  the 
final  image  of  the  point  response  function  for  the  finite-sized  element, 
finite-pa3Sband  array.  See  Fig.  8.  It  is  important  to  emphasize  that  the 
point  response  functions  in  Fig.  8  are  what  a  detector  array  in  our 
"conceptual"  interferometer  would  measure  for  a  point  source. 

While  the  point  response  functions  for  the  type  of  arrays  under  con¬ 
sideration  are  necessarily  two  dimensional,  the  high  degree  of  symmetry  in  the 

L  class  of  arrays  we  studied  means  that  the  average  radial  dependence  of  the 

point  response  function  is  of  interest.  We  used  the  program  POWER  (listed  in 
Appendix  B)  to  read  the  AIPS  map  file  and  average  over  azimuth  to  obtain  the 
profiles  shown  in  Fig.  9.  The  broad  pedestal  in  each  profile  is  the  remnant 
of  the  individual  response. 

1  In  the  case  of  a  complex  source,  the  response  would  be  simply  the  con¬ 

volution  of  the  point  response  function  we  have  calculated  with  the  true 
source.  Deconvolution  to  obtain  the  maximum  possible  resolution  with  high 
dynamic  range  is  possible. 


! 
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Figure  6.  Bandpass  correction  B(x,y)  for  the  15  element  array  considered  in 
detail.  This  corresponds  to  3W.MAP  in  Fig.  4.  Note  that  the  color  wedge  along 
the  bottom  has  been  adjusted  to  exaggerate  the  fainter  parts  of  the  image. 
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Figure  8.  Point  response  functions  for  the  feu-  discrete-el 
apertures  (finite  bandpass  and  finite  elements)  considered 
Upper-left  is  for  a  three  station  interferometer;  upper-rig 
element  array;  loner-left  for  a  15-element  array;  and  lower 
PT-element  array.  Only  one  quadrant  has  been  presented  of  a 
i  function  for  compactness.  The  scale  is  the  same  for  vacs  in 

|  improvement  in  the  point  response  function  as  the  number  of 

1  can  be  seen.  These  maps  correspond  to  PRF.MAP  in  Fig.  U. 


enent  unfilled 
in  detail, 
ht  for  a  seven 
r  Lght  r  >r  i 
ach  point  respons 
age,  so  that  the 
st it  ions  increase 


19 

(Page  20  Blank) 


owe  r 


AVERAGE  POINT  SPREAD  FUNCTIONS 


Figure  9.  Point  response  functions  for  3-,  7-,  15—,  and  27-element  unfilled 
apertures  averaged  in  azimuth  as  a  function  of  radial  distance  from  the 
origin. 


One  approximation  made  in  our  calculations  of  point  response  functions 
was  the  assumption  that  the  individual  elements  were  small  in  comparison  with 
the  overall  size  of  the  array.  In  AIPS,  to  account  for  identical  finite 
elements  exactly  would  require  adjusting  the  convolution  function  to  be  used 
for  gridding  the  data  into  the  UV  plane  so  that  the  gridded  data  would  be  the 
autocorrelation  (weighted  UV  coverage)  of  the  individual  array  elements.  This 
would  be  a  minor  correction  in  the  analysis  we  did. 

RESULTS 

The  point  response  functions  in  Fig.  8  are  the  principal  result  of  this 
report.  The  images  in  Fig.  8  are  the  accurate  two-dimensional  representation 
of  the  point  response  functions  for  the  finite-bandpass,  finite-element, 
unfilled  arrays  which  we  studied  using  the  analysis  of  the  previous  section. 
The  point  response  functions  have  been  shown  using  a  color  wedge  which  sup¬ 
presses  the  central  maximum  (which  appears  only  as  a  few  pixels  of  red  in  the 
lower  left  of  each  image)  in  order  to  make  the  differences  between  the  point 
response  functions  apparent.  The  six  radial  sidelobe  arms  in  the  3-element 
case  and  the  fourteen  radial  sidelobe  arms  in  the  7-element  case  are  remnants 
of  incomplete  sidelobe  suppression  by  the  delay  function  of  the  type  shown  in 
Fig.  6.  As  the  number  of  elements  Increases,  this  delay  sidelobe  supression 
becomes  more  complete. 

In  comparison,  the  point  response  functions  for  the  circular ly- 
symmetric,  monochromatic  cases  we  considered  are  both  simpler  and  much  easier 
to  calculate.  The  two  dimensional  point  response  functions  shown  in  Fig.  8 
are  more  complex,  but  are  easily  understood  in  terms  of  the  circularly  sym¬ 
metric  cases. 

The  sidelobes  of  the  point  response  functions  shown  in  Fig.  8  for  the 
3-.  7-.  15-.  and  27-element  unfilled  apertures  with  finite  bandpass  and  finite 
element  size  have  two  major  components:  (1)  A  pedestal  caused  by  the  response 
of  the  Individual  elements  in  the  array,  and  (2)  Oscillations  about  this 
pedestal  due  to  the  synthesized  beam  component.  The  peak  level  of  the  smooth 
pedestal  component  decreases  as  1/N,  as  can  be  seen  in  Fig.  9  where  the  point 
response  functions  have  been  averaged  in  azimuth.  The  rms  noise  level  of  the 
synthesized  beam  component  also  decreases  as  1/N.  Thus  the  sidelobes  decrease 
rapidly  as  the  number  of  elements  increases  until  approximately  7  to  15  ele¬ 
ments  are  in  the  array.  Larger  numbers  of  elements  then  continue  to  improve 
the  performance,  but  only  at  a  slower  rate;  a  specific  desired  imaging  cap¬ 
ability  may  lead  to  arrays  of  these  dimensions.  At  the  same  time,  the  samp¬ 
ling  in  the  U,V  (Fourier  transform)  space  is  becoming  more  complete  and  the 
synthesized  beam  sidelobes  are  decreasing.  By  the  time  N  is  7,  the  first 
sidelobe  is  19%  but  the  average  oscillations  about  the  1/N  pedestal  level  of 
1 4%  are  at  the  2J  level.  For  the  array  geometries  studied  in  this  report,  the 
larger  number  of  elements  also  decreases  the  minimum  spacing  in  the  U,V  plane, 
lowering  the  lowest  spatial  frequency  sampled,  and  moving  the  first  grating 
sidelobe  (a  circular  ring  for  this  geometry)  further  from  the  center. 
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The  point  response  functions  in  Fig.  8  could  be  improved  through  image 
processing,  with  the  degree  of  improvement  depending  on  the  field  to  be  imaged 
as  well  as  on  the  array  characteristics.  Because  the  array  geometry,  wave¬ 
length,  and  bandpass  would  be  known,  the  point  response  function  could  be  used 
to  deconvolve  the  measured  image  and  yield  an  accurate  representation  of  the 
true  image.  A  limitation  is  that  the  presently  popular  algorithms  for  decon¬ 
volution  yield  better  {higher  contrast)  results  as  the  sidelobes  are  lower, 
and  as  the  measured  signals  fill  a  smaller  percentage  of  the  image  pixels. 

It  is  thus  possible  to  specify  the  required  number  of  elements  for  an 
array  in  two  ways.  First,  the  raw  (no  image-processing)  contrast  ratio  re¬ 
quired  in  the  output  image  could  be  specified.  Since  the  number  of  elements 
determines  the  sidelobe  level  more  or  less  directly  for  this  family  of  un¬ 
filled  arrays,  it  is  possible  to  invert  the  process  and  specify  a  required 
contrast  ratio  for  the  unprocessed  image  and  yield  the  number  of  elements. 

For  instance,  in  a  satellite  surveillance  system  where  image  processing  in 
space  is  prohibited,  a  requirement  of  a  25:1  raw  contrast  specifies  25  ele¬ 
ments.  However,  it  should  be  emphasized  that  most  applications  do  allow  image 
processing,  even  if  not  with  electronic  digital  computers,  and  such  processing 
can  improve  the  images  significantly.  This  suggests  the  second  method  to 
specify  the  required  number  of  elements:  select  the  desired  contrast  ratio  for 
a  processed  image.  Since  the  contrast  in  a  deconvolved,  processed  image 
depends  not  only  on  the  array  configuration  but  also  on  the  complexity  of  the 
image  and  the  signal-strength  in  the  image,  some  broad  characteristics  for  the 
field  being  imaged  would  need  to  be  specified.  Depending  on  the  field  being 
imaged,  10-  to  50-fold  increases  in  contrast  ratio  could  be  achieved.  Thi3 
sort  of  improvement  is  routinely  achieved  in  images  produced  from  radio  inter¬ 
ferometers. 

CONCLUSIONS 

This  report  demonstrates  a  new  technique  for  modeling  point  response 
functions  for  complex  interferometer  systems  taking  into  account  both  finite 
bandpass  and  finite  element  size.  New  interferometer  systems  can  be  fully 
evaluated  with  a  few  weeks  of  effort.  While  the  analysis  system  is  not 
trivial  to  use,  it  i3  straight  forward  as  shown  in  Fig.  4. 

In  view  of  the  limits  to  the  overall  size  of  the  point  response  function 
Imposed  by  the  effects  of  broad  bandwidth  (see  above),  it  is  possible  to 
choose  between  the  two  .basic  types  of  interferometers  used  in  imaging.  The 
first  type,  commonly  used  in  radio- frequency  imaging  arrays  such  as  the  VLA, 
uses  phase  steering  to  map  an  image.  The  maximum  image  size  is  limited  both 
by  the  UV  coverage  of  the  interferometer  (the  usual  case)  and  the  bandpass 
response  of  the  array  and  receivers.  With  the  narrow  bandwidths  available  to 
radio-frequency  detectors,  this  is  a  practical  technique  and  allows  the  map¬ 
ping  of  large  fields  with  high  resolution.  The  second  type  using  delay  steer¬ 
ing  is  required  where  broad-band  phase  stable  correlators  are  impractical  to 
construct,  such  as  at  optical  wavelengths.  The  advantage  of  the  second  method 
(where  beams  from  the  individual  elements  are  combined  coherently  using 
mirrors  and  focused  onto  a  detector  array)  is  that  large  bandwidth  is  easily 
achieved  along  with  a  narrow  point  response  function.  The  detectors  for  such 
an  interferometer  are  directly  measuring  the  convolution  of  the  point  response 
function  with  the  true  image;  for  many  sorts  of  images  reliable  deconvolution 
(restoration)  techniques  exist,  resulting  in  interferometer  performance  which 
equals  that  from  a  filled  aperture. 


The  problem  of  restoring  a  complex  image  is  expected  to  be  worse  for  the 
case  of  a  quasi-optical  system  viewing  the  earth  than  for  the  case  of  the 
radio  images  of  celestial  objects  that  radio  astronomers  typically  restore. 
However,  many  of  the  standard  radio  astronomy  image  processing  techniques  are 
directly  applicable.  It  is  important  to  note  there  are  proven  techniques  to 
deconvolve  images  to  reduce  the  sidelobe  levels  below  those  presented  in  Fig. 
9.  The  exact  effects  of  sidelobe  levels  and  the  extent  to  which  they  could  be 
removed  will  depend  on  the  exact  constraints  and  design  parameters  for  a  real 
system  as  well  as  on  the  characteristics  of  the  field  under  observation. 

As  a  demonstration  of  the  discrete  analysis  system,  a  spaceborne  earth 
monitoring  optical  (350-650  nm)  interferometer  system  was  modeled  to  determine 
the  number  of  elements  needed  to  provide  satisfactory  performance.  We  fixed 
the  maximum  baseline  length  at  10  m,  the  size  of  each  element  at  10  cm,  the 
geometry  to  evenly  spaced  elements  on  a  ring,  and  varied  the  number  of  ele¬ 
ments  between  3  and  99. 

We  find  the  sidelobe  level  for  a  monochromatic  system  varies  as  the 
reciprocal  of  the  number  of  elements  over  the  range  of  interest  in  this  study. 

This  establishes  an  upper  limit  for  any  realistic  system,  because  the  finite 
element  size  reduces  sidelobes  as  was  shown  above. 

We  find  the  realistic  bandpasses  of  interest  limit  the  sidelobes  to  such 
a  level  that  they  are  of  little  concern  for  arrays  with  more  than  a  few  ele¬ 
ments. 


We  find  that  the  difference  between  an  addition  correlator  and  a  multi¬ 
plication  correlator  is  not  significant  in  either  the  difficulty  of  analysis 

or  the  performance  of  a  large  number  of  elements. 
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APPENDIX  A 


BANDWIDTH  EFFECTS 
INTRODUCTION 

The  following  discussion  closely  follows  the  development  of  Thompson 
(1982)  which  was  for  the  particular  case  of  the  VLA  for  radio  astronomy. 

First  the  two  element  interferometer  with  rectangular  bandpass  is  considered. 
Non  rectangular  bandpasses  only  make  the  mathematics  harder  to  solve  in  closed 
form,  but  are  not  intrinsically  difficult.  The  solution  for  two  elements  is  a 
sin(x)/x  function.  In  going  over  to  the  array  case,  a  superposition  of 
solutions  with  the  correct  weighting  is  required.  In  practice  it  is  appro¬ 
priate  to  restart  the  analysis.  However,  this  approach  requires  the 
fractional  bandpass  to  be  small  (^  5% )  in  order  for  the  mathematics  to  reduce 
to  simple  form.  For  wide  bandpasses,  the  mathematics  probably  gets  more 
involved,  and  a  computer  sunraation  is  the  expediant  solution.  The  paper  of 
Thompson  paraphrased  in  part  here  sets  the  framework  for  more  detailed 
analysis  at  a  later  time. 


BANDWIDTH  EFFECT  FOR  2  ELEMENT  INTERFEROMETER 


The  part  of  the  interferometer  in  which  the  signals  are  combined  is  the 
correlator,  which  is  basically  a  voltage  multiplier  and  a  time  averager.  Its 
output  is  the  real  part  of  V’  and  is  related  to  the  input  voltage  waveforms  V1 
and  V2  by: 

R  (V')  =  <  V*(t)  V_  (t+t)>  ( A-1 ) 

e  12 

where  the  angular  brackets  represent  a  time  average  and  we  have  taken  the 
complex  conjugate  of  the  first  waveform. 

Let  I(s,v)  represent  £he  brightness  (intensity)  of  the  target  in  the 
direction  of  unit  vector  s  at  a  frequency  v.  The  signal  pow^r  received  in  a 
bandwidth  dv  from  the  source  dfl  is  A(s)+I(s)  dR  dv,  where  A(s)  is  the 
effective  collecting  area  in  direction  s  which  is  assumed  to  be  the  same  for 
each  interferometer  element. 


If  the  inputs  to  the  correlator  are  represented  by  the  complex  analytic 
signals  the  received  power  is  equal  to  the  square  of  the  modulus,  and  thus: 

Iv^  l2=  I v2I2=  |v  v  !  =  A(s)I(s)dRdv.  (A-2) 

In  the  case  where  the  elements  are  not  identical,  A  is  replaced  by  the  product 
of  the  complex  voltage  responses  of  the  elements.  The  arguments  of  v1  and  v 
are  given  by: 


argCv^  =  2ifvt 
arg[v2]  =  2irv(t+t) 


(A-3) 


where  t  is  the  difference  between  the  geometrical  and  instrumental  delays  in 
Fig.  A-1.  If  B  is  the  baseline  vector  for  the  two  antennas, 
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T  =  T-T .  =  8‘S  C"1  -  T . .  (A-4) 

g  i  1 

Then  from  (A-1)  with  appropriate  integration  we  obtain  the  total  response: 

V,(5,ti)dv  =  r  fs  A(5,v)ei2irv(§*3c”1"T  >dndv  (A-5) 

where  S  indicates  that  the  integral  with  respect  to  Q  is  taken  over  the  sur¬ 
face  of  the  hemisphere.  If  we  take  the  Fourier  transform  of  V*  with  respect 
to  we  obtain: 

V(8,v)  =  /•  V' (§,t. )ei2irvTi  dT 
1 

V(8,v)  =  7-  /°/^VA(s,v)I(s,v)ei2irv®'s/c  dft  dv  (A-6) 

Av  v  -Iv/2 
o 

and  the  integral  over  frequency  depends  upon  the  bandpass*.  Often  over  the 
bandwidth  Av  the  quantities  V,  A  and  I  can  be  assumed  constant  and  we  can 
write  them  as  functions  of  S  or  s  only,  with  the  understanding  that  they  refer 
to  a  center  frequency  vQ.  Thus  (A-6)  simplifies  to 

V(§)  =  /s  A(s)I(s)ei2irvoS*s/c  do  (A-7) 

This  is  the  monochromatic  result  that  is  usually  used  for  narrow  band-pa3ses. 

If  we  assume  A(s,v)  is  a  rectangular  bandpass  from  v  -Av/2  to  v  +Av/2, 
then  the  Fourier  transform  of  a  rectangle  is  a  sin(x)/x  ?unction  yielding 


V(8) 


■L 


sin[irAvS*(  s-s  )/c]  .0  *,+  *  s  . 

o  T/*\  i BlS-S  )/C 

-  A(s)I(s)e  o  o 


d» 


( A-8) 


irAv8,(s— s  )/c 
o 

Here  we  have  assumed  that  the  passband  is  rectangular  with  center  frequency  v 
and  width  Av,  and  s  is  the  direction  for  which  the  instrumental  delay  is  0 
adjusted  to  compensate  for  the  geometrical  delay.  In  the  simple  inter¬ 
ferometer  in  Fig.  A1  the  direction  s  is  indicated  by  the  angle 


9 

measured  with  respect  to  a  plane  normal  to  the  baseline.  Since  B«(s-so) 

B(e-e  )  cos  0  we  see  from  (A-8)  that  for  a  point  source  at  position  0  the 
effec?  of  the°finite  bandwidth  is  to  multiply  the  response  by  the  sine  function 


sin[irAvb(  0-0  ) cos0  c”1] 
o  o 


xAvB( 0-0  )cos0  c 
0  o 


-1 


(A-9) 


•The  sign  convention  for  the  exponential,  here  positive  for  transforming  from 
I  to  V,  varies  between  authors  but  is  inconsequential  in  these  general  discus¬ 
sions. 
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Figure  A1 .  Angles  of  the  incident  radiation. 

As  a  function  of  ®,  the  output  is  a  product  of  the  fringe  pattern  and  the  sine 
function  which  modulates  the  fringes.  The  latter  term  can  be  regarded  as 
resulting  from  the  decorrelation  of  the  signals  produced  by  the  inequality  of 
the  time  delays  for  directions  other  than  ®  .  From  the  Weiner-Khinchine 
theorem  the  autocorrelation  function  of  thesignals  is  the  Fourier  transform 
of  the  power  spectrum,  and  hence  the  rectangular  passband  gives  rise  to  the 
sine  function  term.  This  term  limits  the  response  in  a  manner  that  has  some 
similarity  to  an  antenna  beam,  and  the  result  is  sometimes  referred  to  as  the 
"delay  beam”.  It  is  an  example  of  what  has  been  more  generally  referred  to  as 
space-frequency  synthesis.  Swenson  and  Mathur  (1969)  have  shown  that  there  i3 
an  exact  correspondence  between  frequency  and  antenna  spacing  in  determining 
the  spatial  frequency  response  of  a  pair  of  radio  antennas.  Some  radio  arrays 
have  been  constructed  to  make  particular  use  of  this  principle;  see,  for 
example,  Douglas  et  al.  (1973).  It  is  also  of  particular  importance  in  very 
long  baseline  interferometry  (Robertson  et  al.  1982).  Note  however  thalj  the 
envelope  function  falls  to  half  amplitude  at  (®-®Q)  =  0.6c(B  cos  0oAv)~  .  i.e. 
the  delay  'beamwidth'  depends  upon  the  projected  antenna  spacing  B  cos  eQ. 

For  an  array  with  numerous  element  spacings  the  beamwidth  concept  is 
inappropriate  and  it  is  better  to  examine  the  problem  as  in  the  following 
section. 


BANDWIDTH  EFFECT  FOR  AN  ARRAY 


The  relationship  between  the  measured  visibility,  V(u,v),  and  the  distribu¬ 
tion  of  intensity  of  "brightness,  I(x,y)  can,  with  certain  assumptions,  be  writ¬ 
ten  as 


where 

S(u,v) 

w(u,v) 

b(x,y) 

I(x,y) 

»  * 


V(u,v)S(u,v)w(u,v)  I(x ,y)  *  *  b(x,y)  (A-10) 

sampling  function  that  represents  the  distribution  of  the  data 
points  in  the  (u,v)  plane. 

weighting  function  of  the  data  including  any  applied  tapering. 
S(u,v)w(u,v)  is  the  optical  modulation  transfer  function, 
synthesized  beam. 

brightness  distribution  as  modified  by  the  beams  of  the 
individual  elements, 
denotes  Fourier  transformation, 
denotes  two-dimensional  convolution. 
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The  synthesized  beam  is  defined  as  the  Fourier  transform  of  S(u,v) w(u,v) ,  and 
it  is  equal  to  the  response  to  a  point  source  at  the  (x,y)  origin.  (To  be 
precise,  the  response  to  a  point  source  is  often  defined  as  the  mirror  image 
of  the  beam,  i.e.  b(-x,-y).  In  this  case  S(u,v)w(u,v)  is  entirely  real  so 
b(x,y)  is  symmetrical  and  the  distinction  can  be  ignored.  The  effects  dis¬ 
cussed  here  that  result  in  distortion  of  the  response  to  a  point  source  are 
usually  zero  at  the  (x,y)  origin,  which  is  why  we  have  stated  that  b(x,y)  is 
pthe  response  to  a  point  source  at  that  position.  If  b(x,y)  varies  over  the 
(x,y)  plane  equation  (A-10)  does  not,  of  course,  hold.  In  using  (A-10)  we 
therefore  assume  that  the  various  distorting  effects  to  be  described  are 
negligible;  in  particular  we  shall  apply  (A-10)  only  to  an  incremental 
frequency  bandwidth  dv  centered  at  frequency  v  within  the  receiving  passband. 
Let  (uv,vv)  be  the  spatial  frequency  coordinates  that  correspond  to  frequency 
v.  In  reducing  the  visibility  data  in  the  usual  way  for  continuum  observa¬ 
tions  we  treat  the  data  derived  from  the  full  receiving  passband  as  though  all 
of  the  signal  was  at  the  center  frequency  v  .  Thus,  since  u  and  v  are  mea¬ 
sured  in  wavelengths,  the  assigned  values  f8r  frequency  v  are 


v  u  v  v 
,  0Uv  ov'’ 

'  V  •  V  / 

wavelengths.  For  the  frequency  increment  dv  equation  (A-10)  becomes 


V(- 


V  •  V 


— )S(u,v)w(u,v)^£  (— )2u~,  r^-)*  »  b(x,y)  ( A— 1 1 ) 


V 

o  o 


Here  we  have  used  the  similarity  theorem  of  Bracewell  (1965)  which,  in  one 
dimension,  can  be  stated  as  follows: 


if  g(u)  ^  f(x) ,  then  g(au)^-|4f  f(£) . 

■  3  ■  3 


Note  that  in  ( A— 11)  it  ha3  been  assumed  that  V  and  I  do  not  vary  appreciably 
with  v  over  the  receiving  passband.  Because  we  assign  values  of  u  and  v  that 
do  not  vary  with  v,  it  follows  that  S,  w  and  b  do  not  vary  with  v.  However, 
this  assignment  introduces  a  scaling  error  v/v  in  the  brightness  coordinates. 

The  measured  distribution  is  obtained  by  averaging  the  right-hand  side  of 
(A-11)  over  the  receiving  passband: 


IjCx.y) 


/  P(v)dv 


*  *  b(x,y) 


(A-12) 


where  P(v)  is  the  power  response  of  the  receiving  system.  In  deriving  (A-12) 
we  have  integrated  over  bandwidth  after  Fourier  transformation  to  the  (x,y) 
domain,  but  the  end  result  is  identical  to  integration  over  bandwidth  in  the 
(u,v)  domain  which  is  implicit  in  the  usual  data-processing 

The  variation  of  the  scale  of  x  and  y  with  v  in  (A-12)  is  clearly  a  radial 
effect,  and  averaging  over  the  bandwidth  therefore  results  in  a  radial 
broadening  or  smearing  of  the  detail.  We  see  from  (A-12)  that  the  measured 
brightness  distribution  is  a  modification  of  the  true  brightness  in  which  the 
radial  smearing  is  first  applied  and  then  the  resulting  distribution  is 


» 


» 


> 


» 


» 
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convolved  with  the  synthesized  beam.  The  extent  of  the  smearing  increases  in 
proportion  to  the  distance  from  the  (x,y)  origin.  However,  r inglobe  responses 
of  the  synthesized  beam  are  smeared  to  the  same  extent  as  the  main  beam 
response,  even  if  they  fall  near  the  origin,  as  illustrated  in  Fig.  A2.  Note 


Figure  A2:  Response  to  a  point  source  displaced  from  the  (x,y)  origin  when 
the  synthesized  beam  contains  a  ringlobe.  The  main  beam  and  all 
sidelobes  are  broadened  in  the  radial  direction  by  convolution 
wi'th  the  radially-smeared  point  source. 

that  the  synthesized  beam,  as  we  have  defined  in  above,  does  not  vary  over  the 
map.  However,  the  response  to  a  point  source  does  vary,  because  of  the  radial 
smearing,  and  this  response  is  often  regarded  as  defining  the  beam  for 
practical  purposes.  Statements  about  the  behavior  of  the  'beam'  can  therefore 
lead  to  confusion  if  the  meaning  of  th’e  term  is  not  made  clear. 

The  radial  smearing  reduces  the  amplitude  of  the  response  to  a  source  by 
an  amount  that  depends  upon  the  source  position  and  the  radial  width  of  the 
response.  As  an  example  we  can  calculate  the  reduction  for  a  point  source  for 
a  case  where  the  synthesized  beam  is  represented  by  a  Gaussian  distribution. 
Figure  A3  shows  a  radial  section  through  the  response  to  such  a  source  at  a 


»  4 
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I/0/r,A!/ 


P(l/) 


,'[tr-r,)2/2o’2] 


r1{l-Ai//2i/0)J  '  Lr1(1  +  Ai'/2i/0) 


Figure  A3:  Radial  section  through  point  source  at  r.  which  is  spread  into  a 
rectangular  distribution  by  the  finite  receiving  bandwidth.  The 
Gaussian  function  represents  the  syntehsized  beam. 

distance  r  from  the  origin.  We  shall  represent  the  bandpass  function  P(v)  by 
a  rectangular  function  of  width  Av,  and  ignore  the  factor  (v/v  )*  in  the 
numerator  of  (A-12)  since  it  is  unlikely  that  it  differs  from  unity  by  more 
than  5%.  Then  the  broadened  source  is  represented  by  the  rectangle  in  Fig. 

A3,  and  the  peak  response  by  the  area  under  the  product  of  the  two  function. 
The  ratio,  RQ,  of  the  peak  response  relative  to  that  when  av  tend3  to  zero  is 


-M1 

r.  Av  J 
1  r. 


r1(l+Av/2vQ) 

e  (r  ri)  /2°  dr  =  - , - ^ 

ri(l-Av/2vo)  riav 


erf  ( 


2/Z  a  v 


where  erf  (  )  is  the  normal  error  function.  The  half-power  width  of  the 
synthesized  beam,  eb,  is  2.35a,  and  in  terras  of  eb 


8  v  0 . 832r. Av 

rb  =  L064(n5)erf('vr) 

1  DO 


A  graph  of  this  function  is  shown  in  Fig.  A4.  The  response  is  reduced  by 
a  factor  1/2  when  (r.Av/9bv )  =  2.1,  and  for  r.Av>>9bv  ,  the  response 
decreases  approximately  as  °r.Av)~  .  ' 


E 


Figure  A4:  Decrease  in  the  response  to  a  point  source  with  increasing 
distance  from  the  field  center,  from  equation  (A-14). 


Although  the  peak  response  to  a  point  source  i.e.  the  apparent  brightness 
in  the  map  is  decreased  by  the  radial  smearing,  the  response  is  also  broadened 
and  it  is  easy  to  see  that  the  integrated  flux  density  remains  unchanged.  The 
broadening  of  the  response  for  a  point  source  can  be  calculated  by  convolving 
the  bandpass  function  with  the  synthesized  beam  profile  as  indicated  by 
equation  (A-12) . 
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APPENDIX  B 


PROGRAM  SOURCE  CODE  and  DATA  FILES 


I.  Files  used  in  running  FAKE: 

The  following  DCL  command  procedure  runs  the  Caltech  VLBI  program  FAKE 
for  the  case  of  15  elements  whose  locations  are  specified  in  the 
STATIONS.DAT  file.  The  start  time  has  been  selected  so  as  to  give  one 
baseline  exactly  north-south,  for  compatibility  with  later  programs. 

FAKE. COM 


!  Assignment  to  specify  station  position  file: 

$  ASSIGN  STATIONS.DAT  STATIONS 
$ 

$  FAKE  !  Start  the  program 

OUTPUTS  15STA.DAT'  !  Nane  of  output  data  file 

MODELS' POINT. MOD'  !  Name  of  input  model  file 


YEARs1983 
DAY  356 

START  23:56:00 
STOP  00:03:00 
INTEG  600 


1  Year,  date,  start  and  stop  selected  so 
!  as  to  give  a  single  integration  period 
!  for  the  data  set  with  the  correct 
l  orientation 


FREQ  27037  !  This  frequency,  in  combination  with 

!  the  STATIONS.DAT  file  gives  an  array  of 

!  diameter  20  million  G 

RA  00:00:01  !  Source  position 

DEC  89:50:00  ! 


STATIONS  ' 1 5N 1 '  ,  '15N2',  ' 1 5N3* .  '15N4* ,  '15N5', 
'15N6',  ' 15N7' ,  ' 1 5N  8 '  ,  '15N9',  '15N10', 
'15N11 '  , ' 15N12'  ,'15N13' ,'15N14' ,'15N15' 

PHASES  0.  ERRADD  0.000001  ERRMULT  0.0  / 


Data  files  used  by  FAKE 


The  following  is  an  example  of  the  STATIONS.DAT  file  used  by  FAKE. 
Note  that  in  its  current  incarnation,  FAKE  requires  this  to  be  in  fixed 
format . 

STATIONS.DAT 


5N1 

0. 

89.0 

24.0 

80.0 

I5N2 

0. 

89.0 

48.0 

80.0 

I5N3 

0. 

89.0 

72.0 

80.0 

5N4 

0. 

89.0 

96.0 

80.0 

I5N5 

0. 

89.0 

120.0 

80.0 

5N6 

0. 

89.0 

144.0 

80.0 

I5N7 

0. 

89.0 

168.0 

80.0 

5N8 

0. 

89.0 

192.0 

80.0 

I5N9 

0. 

89.0 

216.0 

80.0 

5N10 

0. 

89.0 

240.0 

80.0 

I5N11 

0. 

89.0 

264.0 

80.0 

I5N12 

0. 

89.0 

288.0 

80.0 

I5N13 

0. 

89.0 

312.0 

80.0 

5N14 

0. 

89.0 

336.0 

80.0 

I5N15 

0. 

89.0 

360.0 

80.0 

The  following  is  an  example  of  the  fake  source  file  POINT. MOD 

POINT. MOD 


!  Sample  fake  data  file 
! 

!  The  following  is  a  1-Jy  point  source. 

!  Flux  Radius  Theta  Axis  Ratio  Phi  Type 
1  0  0  0  0  0  1 


ooo  non  o  o  n 


r 


> 


II.  Bandpass  program 


P 


hi 


C  Program  BAND 

C 

C  BAND  computes  the  bandwidth-dependent  beam  for  an  interferometer 

C  by  summing  the  appropriate  sine  functions  into  an  array.  It 

C  assumes  the  elements  are  evenly  spaced  on  a  ring;  only  odd  numbers 

C  of  elements  up  to  MAXS  elements  are  allowed. 

C 

C  When  linking  this  program,  both  the  Caltech  VLBI  object  library 

C  and  the  Caltech  public  object  library  need  to  be  accessed. 

C 

C  Parameters: 

CHARACTER* ( *)  VERSION 
CHARACTER* 12  ELAPSE 

INTEGER*4  INC , OUTC , OUTDAT, MAXS , MAXB , MAXP , NPARS 
REAL*8  PI.VELC 

PARAMETER  (VERSION=' 1 .0  -  1983  January  28’) 

PARAMETER  ( INC=5 ,0UTC=6 , 0UTDAT=1 0 ) 

PARAMETER  (MAXS=99)  !  maximum  number  of  stations 

PARAMETER  (MAXB=(MAXS-1)»MAXS/2) 

PARAMETER  (MAXP=MAXB)  !  maximum  number  of  u,v-points 
PARAMETER  (NPARS=14)  !  input  parameters 
PARAMETER  (PI=3 . 141 592653589793D0)  !  pi 

PARAMETER  ( VELC=2.997925D8)  !  speed  of  light 

C 

REAL*8  PARS(NPARS) .VALS(NPARS) .ENDMRK 
REAL*4  POSANG(MAXS) ,BUFFER(51 2,512) 

INTEGER*2  IARRAY(512) ,IX, IY 
CHARACTER*64  OUTDSN 
EQUIVALENCE  (OUTDSN, VALS(7)) 


CHARACTER 

T0DAY*9 ,N0W*8 

DATA 

PARS/  ' NSTA' , 

'DIAMETER',  'FOV',  * LMIN' , 

1 

' LMAX' , 

'PIXELS' ,  '0UTFILE' ,  7*' 

DATA 

VALS/  3D0, 

10D0,  2048D0,  3500D0, 

1 

6500D0, 

512D0,  8*'  '  / 

DATA 

ENDMRK/  '/'  / 

►  < 


»  1 


*  i 


>  < 


Introduction 


CALL  DATE(TODAY) 

CALL  TIME (NOW) 

WRITE(0UTC, 1001 )  VERSION 

Control  parameters 

10  WRITE(0UTC, 1001 )  *  Parameters  (type  HELP  for  details):' 
MODE  =  1 

IF  (ISTERM('SYS$INPUT*))  MODE  =  2 

CALL  KEYIN( PARS, VALS, NPARS, ENDMRK, MODE, INC, OUTC) 


Check  parameter  values 


35 


IF  (VALS(1)  .GT.  MAXS  .AND. 

INT(VALS(1 )-1D0)/2  .NE.  INT(VALS(1))  .AND. 

VALS(1)  .NE.  2  .AND.  VALS(1)  .NE.  4)  THEN 

WRITE  (OUTC, 1000)'  Illegal  number  of  stations.  Try  again' , 

' ,  please.' 

GO  TO  10 

ELSEIF  (VALS(6)  .NE.  512  .AND. 

VALS(6)  .NE.  256  .AND.  VALS(6)  .NE.  128  .AND. 

VALS(6)  .NE.  64  .AND.  VALS(6)  .NE.  32)  THEN 

WRITE  (OUTC, 1000) *  Illegal  number  of  pixels.  Try  again,', 

'  please.' 

GO  TO  10 
ELSE 

Attempt  to  open  output  file 
IF  (OUT DSN  .EQ.  '  ' )  GO  TO  20 
CLOSE  (UNIT=10) 

IMSIZE  =  VALS(6) 

OPEN(UNIT=OUTDAT,NAME=OUTDSN,STATUS='NEW' , 

FORM=' UNFORMATTED* ,RECORDTYPE= 'FIXED' ,RECL=IMSIZE/2,ERR=20) 
INQUIRE(UNIT=OUTDAT, NAME=OUTDSN) 

GOTO  30 

WRITE  (OUTC, 1000)'  Unable  to  open  output  file.  Please  try', 
'  again.' 

GO  TO  10 
CONTINUE 
ENDIF 

NSTA  =  VALSO) 

Write  out  input  parameters  with  units 
WRITE  (OUTC, 1002)  NSTA, (VALS(I) ,1=2,5) , IMSIZE, 

OUTDSN(  1 :  LEN1  (OUTDSN) ) 

Calculate  array,  one  position  angle  and  baseline  at  a  time 
IF  (NSTA  .EQ.  2)  THEN 
POSANG(I)  =  90 
ELSE 

DO  1=1 ,(NSTA-1)/2 

POSANG ( I )  =  I* (180. /NSTA)  -  90. 

ENDDO 

ENDIF 

NBASE  =  NSTA* (NSTA-1 ) /2 
N LENGTH  =  ( NSTA-1 )/2 

Begin  array  calculations  here 

DELAY  =  1D10  *  PI  »  ABS(1/VALS(4)  -  1/VALS(5)) 

FOV  =  (VALS(3)  *  PI)  /  (180.  *  3600.  *  1000) 

CALL  RCLOCK (ELAPSE, CPU) 

WRITE  (6,1000)  '  Ready  to  start  long  calculation.' 

WRITE  (6,1003)  ELAPSE, CPU 


o  o  non  n  no  on 


DO  IL  =  1 , NLENGTH 

BSLINE  =  VALS(2)*C0S( (IL-0.5) *PI/NSTA) 

DO  IB  =  1 ,(NSTA-1 )/2 

C  Calculate  1st  half  of  array: 

CALL  ARRAY(DELAY, BSLINE, FOV, POSANG(IB) .IMSIZE, BUFFER) 

ENDDO 

ENDDO 

C 

CALL  RCLOCKCELAPSE, CPU) 

WRITE  (6,1000)  '  Long  calculation  finished.' 

WRITE  (6,1003)  ELAPSE, CPU 
C 

C  Now,  load  up  upper  half  of  array 

DO  IX  =  1,  IMSIZE 

DO  IY  =  2,(IMSIZE/2)+1 

BUFFER ( IMS IZE-IX, IMSIZE+2-IY)  =  BUFFER(IX.IY) 

ENDDO 
ENDDO 

Multiply  this  first  half  of  array  by  2  for  normalization  reasons 
DO  IY  =  1,  IMSIZE 
DO  IX  =  1,  IMSIZE 

BUFFER(IX.IY)  =  BUFFER(IX, IY)  +  BUFFER(IX, IY) 

ENDDO 
ENDDO 

Add  in  the  pa=90  and  pa=0  contribution  (if  any) 

DO  IL  =  1, NLENGTH 

BSLINE  =  VALS(2)#COS( (IL-0.5) *PI/NSTA) 

CALL  ARRAY (DELAY, BSLINE, FOV, 90. , IMSIZE, BUFFER) 

IF  (NSTA  .EQ.  (NSTA/2)*2)  THEN 
Nunber  of  stations  is  even 

CALL  ARRAY(DELAY, BSLINE, FOV, 0. , IMSIZE, BUFFER) 

END  IF 
ENDDO 

Now,  flip  array  over  and  add  into  final  result,  then 
write  out  array  into  file  after  normalizing 
PEAK  =  2  *  NBASE 
DO  IY  =  1,  IMSIZE 
DO  IX  =  1  .IMSIZE- 1 
IARRAYdX)  =  32767  * 

(BUFFER(IX, IY)+BUFFER(IMSIZE-IX, IY))  /  PEAK 

ENDDO 

IARRAYdMSIZE)  =  0  !  Since  flip  was  not  valid  for  that  column 
WRITE  (10)  (IARRAYdX),  IX  =  1,  IMSIZE) 

ENDDO 

CLOSE  (UNIT=10) 

STOP 


1000  FORMAT(A) 

1001  FORMATC  BAND  calculates  the  delay  beam  for  simple  arrays'  , 

+  /'  (Version  ',A,')'/) 
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1002  FORMAT  ('  Input  parameters: '/ 

+  *  ',I8,2X,'s  nunber  of  stations.*/ 

+  *  ',F8.2,2X,'s  array  diameter  in  meters.’/ 

+  *  ' ,F8.2,2X,*s  field  of  view,  in  mas.*/ 

+  '  ',F8.1,2X,'=  lower  wavelength,  angstroms.*/ 

+  *  ',F8.1,2X,'=  upper  wavelength,  angstroms.'/ 

+  *  '  ,I8,2X,*s  output  image  size.'// 

+  '  Output  to  file  ' ,A) 

1003  FORMAT  ('  Elapsed  time:  '.A,'  CPU:',F15.3> 

C 

END 


C 

C 

C 

C 


C 

C 


a 


c 


> 


► 
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SUBROUTINE  ARRAY(DELAY, BSLINE, FOV, POSANG, IMSIZE, BUFFER) 

Calculate  array:  Delayspi*  d  nu  *  /  c,  in  meters 
bslinesbaseline  in  meters 
fovsfield  of  view,  in  radians 
lmsize=array  size,  in  pixels 

REAL*8  PI 

PARAMETER  (PI=3 . 141 592 653 5897 93D0)  1  pi 

REAL»4  ROW(512) .COLUMN (51 2) .BUFFER (51 2,512) 

INTEGERS  IX,  IY,  ICX,  ICY 

Center  for  AIPS 

ICX  s  IMS  n E/2 

ICY  s  ICX  ♦  1 

R POSANG  =  P0SANG*PI/180. 

XUNIT  s  DELAY  •  BSLINE  «  FOV  /  IMSIZE 


Calculate  array,  with  special  cases  for  0  and  90  degrees 
IF  (POSANG  .EQ.  90.)  THEN 
DO  IYsl ,  IMSIZE 

First,  calculate  a  column  in  the  array  (x=1) 

X  s  XUNIT  *  RADIUS  •  ABS(ICY  -  IY) 

IF  (X  .EQ.  0.)  THEN 
SINCX  =  1 
ELSE 

SINCX  =  SIN(X)  /  X 
ENDIF 

COLUMN(IY)  =  SINCX 
ENDDO 

DO  IYsl,  IMSIZE 
DO  IXsI.IMSnE 

BUFFER(IX, IY)  s  BUFFERCIX, IY)  +  COLUMN(IY) 

ENDDO 

ENDDO 

ELSEIF  (POSANG  .EQ.  0.)  THEN 
DO  IX  s  1 ,  IMSIZE 

First,  calculate  first  ROW  in  array  (y=1) 

X  s  XUNIT  »  RADIUS  *  ABSUCX  -  IX) 


4 


4 


< 
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IF  (X  .EQ.  0.)  THEN 
SINCX  =  1 
ELSE 

SINCX  =  SINCX)  /  X 
ENDIF 

RCW(IX)  =  SINCX 
ENDDO 

DO  IY=1 ,  IMSIZE 
DO  IX=1,  IMSIZE 

BUFFER(IX, IY>  =  BUFFERCIX, IY)  +  RCW(IX) 

ENDDO 

ENDDO 

ELSE 

Calculate  lower  half  of  array 
DO  IX  =  1,  IMSIZE 

First,  calculate  first  row  in  array  (y=1) 

IDX  =  IX-ICX 
IDY  =  ICY-1 

RADIUS  a  SQRT(FLOAT(IDX**2  +  IDY**2)) 

ANGLE  a  RPOSANG  -  ATAN2 (FLOAT (IDX) .FLOAT (IDY) ) 

X  =  ABS(SIN( ANGLE)  »  XUNIT  «  RADIUS) 

IF  (X  .EQ.  0.)  THEN 
SINCX  =  1 
ELSE 

SINCX  =  SIN(X)  /  X 
ENDIF 

BUFFERCIX, 1)  a  BUFFERCIX, 1)  +  SINCX 

DO  IY  a  2. ICY-1 
IDY  a  IY-ICY 

RADIUS  a  SQRT(FL0AT(IDX»»2  +  IDY»*2)) 

ANGLE  a  RPOSANG  -  ATAN2 (FLOAT (IDX) .FLOAT (IDY) ) 
X  a  ABS( SIN (ANGLE)  •  XUNIT  *  RADIUS) 

IF  (X  .EQ.  0.)  THEN 
SINCX  a  1 
ELSE 

SINCX  a  SIN(X)  /  X 
ENDIF 

BUFFER(IX, IY)  a  BUFFERCIX, IY)  +  SINCX 
ENDDO 
ENDDO 


DO  IX=1 , ICX-1 
RADI US= ICX-IX 
ANGLE  a  RPOSANG  -  PI/2. 

X  a  ABS(SIN( ANGLE)  *  XUNIT  *  RADIUS) 

IF  (X  .EQ.  0.)  THEN 
SINCX  a  1 
ELSE 

SINCX  a  SIN(X)  /  X 
ENDIF 

BUFFERCIX, ICY)  a  BUFFERCIX, ICY)  +  SINCX 


O  O  OOO  DO  o  o  ooooooo 


III.  Progran  to  find  radial  point  response  function  profiles. 
Program  POWER 

This  progran  finds  the  power  for  a  512  on  a  side 
AIPS  map,  with  the  output  in  a  1-d  array  with  average  power  as 
a  fn  of  radius  from  256,257  in  the  output  file.  The  output  file 
is  in  LISTING  format  so  that  it  may  be  printed,  etc. 

CHARACTERS  OUTDSN,  INDSN 
PARAMETER  IMSEE=512 
REAL*4  POWER  (IM  SIZE) 

INTEGER*2  INDATA(  IMSIZE,  IMSIZE) , MIN,  MAX,  IX,  IY,  IREC,  IB  IN, 

+  WEIGHT  (IMS  EE) 

Attempt  to  open  files 

10  WRITE  (6,1000)  '$Filename  of  AIPS  map:  ' 

READ  (5, 1000,END=99,ERR=10)  INDSN 
CLOSE  (UNITslO) 

0PEN(UNIT=1 0, NAME=INDSN, STATUS= 'OLD' , READONLY, 

1  FORM=' UNFORMATTED' ,RECORDTYPE= 'FIXED' ,RECL=128) 
INQUIRE(UNIT=10,NAME=INDSN) 

GO  TO  30 

20  WRITE  (OUTC, 1000) '  Unable  to  open  input  file.  ' 

GO  TO  10 
30  CONTINUE 

Get  output  file 

40  WRITE  (6,1000)  ’$Filenane  of  output  file:  ’ 

READ  (5,1000,END=99,ERR=40)  OUTDSN 
CLOSE  (UNIT=1 1 ) 

OPEN(UNIT=1 1 , NAME=OUTDSN, STATUS= 'NEW'  , 

1  F0RM= 'FORMATTED' , CARR IAGEC0NTR0L= 'LIST' , ERR =50) 

INQUIRE(UNIT=1 1 ,  NAME=OUTDSN) 

GO  TO  60 

50  WRITE  (OUTC, 1000)'  Unable  to  open  output  file.  ' 

GO  TO  40 
60  CONTINUE 

WRITE  (6, 1001)'0The  input  file  is  '  ,INDSN(1 :LEN1( INDSN)) 

WRITE  (6,1001 ) 'OThe  output  file  is  ' .OUTDSN ( 1 :LEN1 (OUTDSN) ) 

READ  IN  DATA 
DO  IY  =  1 ,  IMS  EE 

DO  IREC=1 ,  IMSEE/256 

READ  (10)  (INDATA(IX.IY) ,  IX=(IREC-1 ) *256*1 , IREC*256) 

ENDDO 

ENDDO 

Find  max  and  min  of  array 

MIN  =  32767 

MAX  =  INDATA(256,257) 

DO  IY  =  1 ,  IMS  EE 
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DO  IX  =  1 ,  IMSIZE 

IF  (INDATA(IX.IY)  .LT.  MIN)  THEN 
MIN  =  INDATA(  IX,  IY) 

ELSE  IF  (INDATA(IX.IY)  .GT.  MAX)  THEN 
MAX  =  INDATA(  IX, IY) 

ENDIF 

ENDDO 

ENDDO 

MAXI  4  =  MAX 
MINI4  =  MIN 
C 

C  Find  radii  and  sun  up,  with  peak  normalized  to  1 

ICX  =  IMS  IZ  E/2 
ICY  =  ICX+1 

SCALE  =  1  ,/FLOAT(MAXI4  -  MINI4) 

DO  IY  =  1,  IMSIZE 
DO  IX  =  1,  IMSIZE 

RADIUS  =  SQRT(FLOAT(  (IX-ICX)»»2  +  (IY-ICY)««2)) 

IB IN  =  RADIUS  +  1 

POWER (IB IN)  =  POWER (IB IN)  +  (INDATA( IX, IY)-MINI4)«SCALE 
IWEIGHT(IBIN)  =  IWEIGHT(IBIN)  +  1 
ENDDO 
ENDDO 


Write  into  output  file  a  single  vector 
WRITE  (11,1002)  ( POWER  (IX)  .WEIGHT  (IX)  ,  IX  =  1,  IMSIZE) 
C 

CLOSE  (UNIT=1 1 ) 

CLOSE  (UNIT=10) 

CLOSE  (UNIT=12) 

STOP 
C 

99  STOP  'Error:  end  of  input' 

C 

1000  FORMAT(A) 

1001  F0RMAT( A, A) 

1002  FORMAT  (512(F10.6, 110/) ) 

END 
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IV.  Mirror  Pogram 

This  pogram  finds  the  beam  for  an  aperture  of  specified  size  at 
a  given  wavlength. 


PRO  MIRROR, INMAP, OUTMAP, DIAM, FOV 

;  This  procedure  corrects  the  beam  pattern  (INMAP)  for  the  finite 

;  size  in  cm  (DIAM)  of  the  elemental  mirrors. 

;  The  input  (and  output)  maps  have  a  field  of  view  =  FOV 

• 

;  Find  the  array  size  using  the  input  array 

NS  =  SIZE( INMAP) 

IF  (NS(O)  LT  2)  THEN  GOTO, DIM ERR 
IF  ( NS ( 1 )  NE  NS(2) )  THEN  GOTO, DIM ERR 
MAX  =  NS(1 ) 

CENT  =  MAX/2  -  1 
5 

CX  =  CENT 
CY  =  CENT  +  1 

OUTMAP  =  FLTARR(NS(1) ,NS(2)) 

t 

;  The  wavelength  is  fixed  at  500  nanometers  = 

C  =  1.25664E5  *  DIAM 

t 

;  The  constant  SI  is  the  FOV  (mas)  expressed  in  pixels 
;  4.8481 4E-9  =  4.8E-6«E-3  for  mas 

SI  =  4.84814E-9  «  FOV  /  2. 

SI  =  SI  /  MAX 
CSI  =  C  »  SI 


;If  not  2-D  array 
;flag  a  dimension  error 

;This  is  the  center 
because  IDL  counts  from  zero 


;Creat  the  output  array 

0. 5umeter 

;DIAM  is  in  cm 


The  output  array  is  filled  in  one  part  and  copied  to  reduce  the 
number  of  page  faults  on  the  VAX  computer  over  what  we  would  have 
if  the  output  array  were  filled  directly. 

FOR  J  =  1 ,CY  DO  BEGIN 
FOR  I  =  0,J  DO  BEGIN 


;  First  get  the  radial  distance  from  the  center  =  R 
R  =  SQRT(FLOAT( (I-CX) ~2  +  (J-CY)"2)) 

OUTMAP(I.J)  =  BESJKR  *  CSI)  /  (R*SI)  {Calculate  the 

first  order  Bessel  function 
OUTMAP(MAX-(I+2) ,J)  =  0UTMAP(I,J)  ;Use  symmetry  to  store 

END 


END 

FOR  J  r 
FOR  I 


END 

END 

FOR  J  = 
FOR  I 


1 , CY  DO  BEGIN 
=  0,J  DO  BEGIN 

0UTMAP(MAX-(I+2),MAX>-J)  =  OUTMAP(I.J) 
OUTMAP ( I, MAX -J)  =  0UTMAP(I,J) 


1  , CY  DO  BEGIN 
=  0,J  DO  BEGIN 

0UTMAP(J -1,1+1)  =  0UTMAP(I,  J) 

0U1WAP( MAX- ( J+1 ) , 1+1 ) rOUTMAP ( I, J ) 
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RETURN 

DIMERR: 

END 


END 

END 

FOR  J  =  1,CY  DO  BEGIN 
FOR  I  =  O.J  DO  BEGIN 

OUTMAP(MAX-(J+1 ) ,MAX-(I+1 ) )=OUTMAP(I, J) 
OUTMAP ( J - 1 , MAX- ( 1+1 ) ) =OUTMAP ( I , J ) 


END 

END 

NOW  WE  MUST  DO  THE  4  SIDES  FOR  J=0 
FOR  I  =  O.CX  DO  BEGIN 

R  =  SQRT(  FLQAT(  ( I-CX)  *2  +  CY»CY)) 
OUTMAP (1,0)  =  BESJ1 (R#CSI)  /  (R*SI) 
0UTMAP(MAX-(I+2) ,0)  =  OUTMAP(I.O) 
0UTMAP(MAX-1 ,1+1 )  =  OUTMAP(I.O) 

0UTMAP(MAX-1 ,MAX-(I+1 ))=0UTMAP(I,0) 

END 


OUTMAP ( CX, CY)  =  OUTMAPC  CX+1 , CY) 

OUTMAP  =  OUTMAP  •  OUTMAP 
OUTMAP  =  OUTMAP  /  OUTMAP(CX, CY) 

OUTMAP  =  OUTMAP  •  IMAP  {Multiply  the  input  map 


»  i 


»  i 


PRINT,' Input  map  is  not  square  and  at  least  2-D' 
RETURN 


» 


» 


> 


» 


» 


44 


> 


